AD-A147  258  HEAT  TRANSFER  DUE  TO  FILM  CONDENSATION  ON  VERTICAL  i/1 

FLUTED  TUBES(U)  NAVAL  POSTGRADUATE  SCHOOL  MONTEREV  CA 
V  K  GARG  ET  AL.  JUL  84  NPS69-84-006 


UNCLASSIFIED 


F/G  13/1 


NL 


AD- A 147  258 


NPS69-84-006 


NAVAL  POSTGRADUATE  SCHOOL 

Monterey,  California 


Q_ 

O 


!_i_J 
- 1 

€  a 

P 

r:"' 


HEAT  TRANSFER  DUE  TO  FILM  CONDENSATION 
ON  VERTICAL  FLUTED  TUBES 

V.  K.  Garg 
P.  J.  Marto 
July  1984 


Final  Report  for  Period 
January  1983  -  September  1983 


Approved  for  Public  Release;  Distribution  Unlimited 

Prepared  for: 

Naval  Postgraduate  School 
Monterey,  California  93943 


84  11  01  109 


NAVAL  POSTGRADUATE  SCHOOL 
Monterey,  California 


t 

D.  A.  Schrady 
Provost 


Commodore  R.  H.  Shumaker 
Superintendent 


The  work  reported  herein  was  supported  in  part  by  the  Foundation  Research 
Program  of  the  Naval  Postgraduate  School  with  funds  provided  by  the  Chief  of 
Naval  Research. 

Reproduction  of  all  or  part  of  this  report  is  authorized. 

This  report  was  prepared  by: 


Vi jay  K.  Garg 
Adjunct  Professor  of 
ffechanical  Engineering 


P.  J.  Mart(o 

Chairman,  Department  of 
Mechanical  Engineering 


'I 


Reviewed  by: 

Paul  J.  Marto 
Chairman,  Department  of 
Mechanical  Engineering 


Released  by: 


John  N.  Dyer 

Dean  of  Science  and  Engineering 


SECURITY  CLASSIFICATION  OF  THIS  HACK  r*H«  Oali  tnlmfd) 


REPORT  DOCUMENTATION  PAGE 


NS FONT NUMBER 

NPS69-84-006 


4.  TITLE  (mtd  Mill It) 

Heat  Transfer  Due  to  Film 
Condensation  cn  \fertical  Fluted  Hibes 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


ClFiCM'S  C AT  ALOO  NUMBER 


S.  TyFS  OF  FIFO  NT  •  PERIOD  COVERED 


Final  Report 


#.  PERFORMING  oro.  refort  NUMBER 


T.  AUTHOFf«> 


V.  K.  Garg  and  P.  J.  Marto 


S.  FERFORMINS  ORGANIZATION  NAME  ANO  AOORESS 

Naval  Postgraduate  School 
Monterey,  CA  93943 


It.  CONTROLLING  OFFICE  NAME  AnO  ADDRESS 


MONITORING  AGENCY  NAME  A  AOORESVif  SHWMI  IfFM  CmHnlllnt  OfflcFj  IS.  SECURITY  CLASS,  f mt  Ala  '•+•*) 

Unclassified 

HZ'  IRC assi fixation/ downgrading 
scneoule 


IS.  DISTRIBUTION  STATEMENT  |ol  Ala  RapMU 


Approved  for  public  release;  distribution  unlimited 


17.  DISTRIBUTION  STATEMENT  (mi  thm  mkmttmct  mntmrmdln  fl«c*  79,  It  Mtttmrmi  trmm  Rmmmri) 


If.  KEY  WOROS  (Conttnmm  on  rmwmrmm  mUtm  it  nnmmmmn ry  m%4  Itmnttty  kf  ktmck  mmtktr) 

Film  Condensation,  Fluted  Tubes,  Surface  Tension  Forces ,  Heat  Transfer 


20  ABSTRACT  (Continue  an  raaaraa  alVa  II  nacaaaao’  an4  Idnniltr  *F  WaaA  "MlnJ 

An  analysis  of  film  condensation  on  a  vertical  fluted  tube  has  been  made 
considering  gravitational  and  surface  tension  effects  over  the  entire  fluted 
surface,  and  using  surface-oriented  coordinates.  For  the  first  time  surface 
tension  effects  are  determined,  as  they  should,  from  the  shape  of  the 
condensate-vapor  interface  rather  than  the  shape  of  the  flute. 

Two-d imens ional  conduction  within  the  condensate  film  as  vrell  as  in  the  fluted 
tube  is  considered.  A  finite-difference  solution  of  the  highly  non-linear 
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20.  partial  differential  equation  for  foe  film  thickness  is  coupled 

with  a  finite-element  solution  of  the  conduction  problem.  The  procedure 
has  been  tested  on  a  sinusoidal  flute  with  amplitude  to  pitch  ratio  ~0.2. 
A  linear  extrapolation,  cn  a  log-log  basis,  of  our  results  shows  good 
comparison  with  experimental  data. 
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SUMMARY 


An  analysis  of  film  condensation  on  a  vertical  fluted  tube  has  been  made 
considering  gravitational  and  surface  tension  effects  over  the  entire  fluted 
surface,  and  vising  surface-oriented  coordinates.  For  the  first  time  surface 
tension  effects  are  determined,  as  they  should,  from  the  shape  of  the 
condensate-vapor  interface  rather  than  the  shape  of  the  flute. 

Two-dimensional  conduction  within  the  condensate  film  as  well  as  in  the  fluted 
tube  is  considered.  A  finite-difference  solution  of  the  highly  non-linear 
partial  different  _al  equation  for  the  film  thickness  is  coupled  with  a 
finite-element  solution  of  the  conduction  problem.  The  procedure  has  been 
tested  on  a  sinusoidal  flute  with  amplitude  to  pitch  ratio ^.2.  A  linear 
extrapolation,  on  a  log-log  basis,  of  our  results  shows  good  comparison  with 
experimental  data,. 
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1.  INTRODUCTION 


The  U.S.  Navy  has  a  continued  interest  in  the  reduction  of  the  size  and 
I  weight  of  propulsion  components  aboard  both  surface  vessels  and  submarines. 

Some  studies  are  currently  underway  in  the  Navy  to  determine  the  savings  in 
both  weight  and  volume  which  can  occur  if  the  condenser  is  designed  to  use 
1  corrugated  or  indented  tubing  to  enhance  heat  transfer  on  both  the  shell  and 

tube  sides.  While  a  comprehensive  research  program  has  been  underway  at  the 
Naval  Postgraduate  School  to  study  various  enhancement  techniques  for 
horizontal  condenser  tubes  (l,2]  ,  am  attractive  alternative  aboard  submarines 
is  a  vertical  condenser  oriented  aft  of  the  steam  turbine  as  shown  in  Figure 
1(a).  With  this  arrangement,  there  would  be  a  shorter  machinery  stack  length 
compared  to  a  horizontal  condenser  mounted  underneath  the  turbine.  This  would 
allow  a  smaller  diameter  submarine  to  be  built  at  a  significant  reduction  in 
cost.  The  ultimate  aim  of  this  project  is  therefore  to  assess  the  reduction 
in  condenser  weight  and  size  that  may  be  feasible  if  the  vertical  tubes  in 
this  condenser  are  fluted  on  the  outside  (i.e.,  the  steam  side).  Figure  1(b). 

Many  methods  (3l  for  enhancing  condensation  heat  transfer  have  been 
proposed.  Among  them,  Gregorig  [U]  first  recognized  the  importance  of  surface 
tension  in  film  condensation  on  vertical  fluted  tubes.  Thereafter,  many 
experimental  studies  (5-ll)  on  vertical  fluted  surfaces  were  made  to  confirm 
his  findings.  Lustenader,  et  al.,  (5!  employed  a  fluted  vertical  tube  and 
obtained  about  four  times  larger  heat  transfer  coefficients  than  those  on  a 
vertical  smooth  tube.  Camavos  (61  carried  out  experiments  on  a  doubly  fluted 
vertical  tube  and  found  the  augmentation  ratio  of  condensation  heat  transfer 
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rate  to  be  about  six.  Thomas  (7,8]  tested  a  vertical  tube  with  longitudinal 
rectangular  fins  and  wires.  The  augmentation  ratio  was  about  eight.  Newson 
and  Hodgson  19 1  manufactured  various  condenser  tubes  of  highly  enhanced  heat 
transfer  performance.  Combs  and  his  co-workers  (10,  11]  at  Oak  Ridge  National 
Laboratory  found  the  augmentation  ratio  to  be  around  five  for  condensation  of 
ammonia  and  refrigerants  oa  vertical  fluted  tubes. 

Gregorig's  theoretical  model  has  also  been  improved  upon  but  much  still 
remains  to  be  done.  Edwards,  et  al . ,  [12]  proposed  a  condensation  model  on  a 
heat  transfer  surface  of  triangular  fins  on  the  assumption  of  liquid  film 
attaching  to  the  tip  of  fin  with  a  contact  angle.  The  effect  of  a  locally 
thin  condensate  film  on  the  side  of  the  fin  (Figure  2)  was  not  considered  by 
them  as  well  as  by  Fujii  and  Honda  [13] •  The  latter,  however,  did  consider 
two-dimensional  conduction  within  the  condensate  film  in  the  trough  region 
(Figure  2)  and  within  the  fin  wall.  Mori,  et  al.,  IlU]  considered  the  effect 
of  a  thin  condensate  film  on  the  side  of  the  fin  but  neglected  the  variation 
of  its  thickness  in  the  vertical  direction.  They  also  considered  only 
one-dimensional  conduction  within  the  fin  and  neglected  any  conduction  within 
the  condensate  film  in  the  trough  region.  Hirasawa,  et  al.,  [ 15 1  improved 
upon  [14]  in  that  they  included  the  variation  of  a  thin  condensate  film  in  the 
vertical  direction  but  neglected  conduction  within  the  fin  and  the  film 
completely.  Panchal  and  Bell  (l6,  17]  also  neglected  conduction  within  the 
fin  and  the  film  while  analysing  a  sinusoidal  fluted  tube,  but  later  found 
that  two-dimensional  conduction  is  important  within  the  fin  and  the  film  for  a 
triangular  fin  (l8] .  A  recent  empirical  analysis  by  Barnes  and  Rohsenow  |l9l , 
based  largely  on  (l6,  20],  reports  an  augmentation  ratio  of  about  fifteen  for 
condensation  of  steam  on  a  fluted  surface  while  most  earlier  studies 


considered  condensation  of  a  refrigerant  and  found  much  smaller  augmentation 
ratios.  With  present  day  knowledge,  it  is  difficult  to  ascertain  whether  this 
discrepancy  is  due  to  different  fluids  or  due  to  a  questionable  analysis. 

All  theoretical  analyses  discussed  above  break  up  the  fluted  surface  into 
basically  two  parts.  In  the  portion  near  the  crest,  gravity  is  neglected  in 
comparison  to  the  surface  tension  effect,  while  in  the  portion  near  the 
trough,  the  reverse  is  done.  The  patching  between  the  two  regions  is  carried 
out  at  a  point  that  is  selected  quite  arbitrarily  at  times  [15—18 ] .  This 
isolation  of  the  two  important  effects,  gravity  and  surface  tension,  is 
Justified  on  the  basis  that  the  condensate  film  is  thick  in  the  trough  region 
and  thin  over  the  crest.  This  is,  however,  not  true  over  the  initial  portion 
of  the  tube  length  in  the  vertical  direction  nor  during  the  initial  portion  of 
the  tube  Just  below  a  condensate  drainage  skirt  (Figure  1(b)).  These  initial 
portions  may  be  a  few  pitches  in  length.  A  recent  analysis  by  Stack  and 
Merkle  (21 ]  does  attempt  to  solve  the  complete  equation  with  both 
gravitational  and  surface  tension  effects  included  over  the  entire  flute  but 
has  three  major  drawbacks.  First,  when  the  condensate  film  in  the  trough 
region  has  thickened,  the  analysis  treats  it  in  the  boundary-layer-sense,  that 
is,  it  neglects  all  velocity  gradients  except  those  in  the  direction  normal  to 
the  fluted  surface.  Second,  the  analysis  is  restricted  to  impracticably  low 
values  of  amplitude  to  pitch  ratio  of  the  flute  (0.02  and  0.0k)  owing  to  the 
use  of  a  Cartesian  coordinate  system  rather  than  a  surface-oriented  coordinate 
system.  Third,  the  analysis  does  not  consider  any  heat  transfer  effects. 

Their  analysis  is  only  confined  to  finding  the  condensate-vapor  interface. 
Fujii  and  Honda  [ 13 ]  also  considered  the  entire  thin  film  as  one  piece  over  an 
initial  length  of  the  tube  (about  1/8  of  the  pitch)  but  neglected  the 
important  surface  tension  effect. 


A  major  deficiency  of  all  the  above  theoretical  analyses  except  I 21 1  lies 
in  the  way  surface  tension  effect  is  determined.  Contrary  to  the  general 
belief,  this  effect  does  not  depend  upon  the  curvature  of  the  condensate-vapor 
interface.  Instead  it  depends  upon  the  variation  of  this  curvature  along  the 
interface.  However,  since  the  location  of  this  interface  is  unknown  a  priori, 
many  analyses  (12,  13 1  simply  determine  it  on  the  basis  of  the  known  flute 
shape  owing  to  the  argument  that  in  the  crest  region  where  the  surface  tension 
effect  is  important,  the  film  is  very  thin.  Such  analyses  cannot  be  applied 
at  all  to  triangular  or  rectangular  fins  since  their  curvature  as  well  as  the 
variation  of  curvature  along  the  flute  is  zero.  The  same  is,  however,  not 
true  for  the  condensate-vapor  interface.  Even  for  sinusoidal  flutes,  there 
are  large  differences  between  the  curvature  and  its  variation  along  the  curved 
surface  for  the  given  flute  and  the  actual  interface  (see  Fig.  10  and  its 
discussion  in  Section  5).  Moreover  those  analyses  that  do  claim  to  determine 
the  surface  tension  effect  based  on  the  actual  shape  of  the  condensate-vapor 
interface  are  also  questionable  in  their  claim  as  shown  in  Appendix  A. 

Another  problem  with  all  earlier  theoretical  analyses  except  (13)  is  that 
conduction  within  the  condensate  film  and  the  tube  wall  is  either  completely 
neglected  or  considered  to  be  at  most  one-dimensional.  Panchal  and  Bell  [l8] 
point  out  clearly  that  two-dimensional  conduction  should  be  considered  within 
both  the  condensate  film  and  the  tube  wall. 

In  an  attempt  to  take  care  of  these  deficiencies,  an  analysis  has  been 
developed  during  this  study  that  solves  the  complete  equation  with 
gravitational  as  well  as  surface  tension  effects  included  over  the  entire 
fluted  surface  using  the  surface-oriented  coordinate  system.  In  addition,  for 
the  first  time,  surface  tension  effects  are  determined  as  they  should,  from 


the  actual  shape  of  the  condensate-vapor  interface  rather  than  the  shape  of 
the  flute.  Two-dimensional  conduction  within  the  condensate  film  as  well  as 
the  fluted  tube  is  also  included.  A  finite-difference  solution  of  the  highly 
non-linear  partial  differential  equation  for  the  condensate  film  thickness  is 
coupled  with  a  finite-element  solution  of  the  conduction  problem.  Details  of 
the  analysis  follow  in  subsequent  sections. 
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2.  ANALYSIS 


Considering  the  condensate  as  a  viscous,  incompressible  Newtonian  fluid, 
and  setting  up  a  curvilinear  orthogonal  coordinate  system  (x^,  x^,  z)  as  shown 
in  Fig.  3,  the  continuity  equation  yields  [22] 


F.  +  x^  3x^ 


2 

R  +  x. 


(2.1) 


where  (uj  ,  u2,  w)  are  the  velocity  components  in  the  (x  ,  x2,  z)  directions 
respectively,  and  R(x1)  is  the  radius  of  curvature  of  the  fluted  tube.  The 
momentum  equations  are  (cf.  Ref.  22,  p.  68) 
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or  3/3z,  and  u2  «  or  w.  This  simplifies  the  momentum  equations  (2.2)  - 
(2.U)  to 


1  R  3p 
U  R  +  x2  3x1  ’ 


(2.5) 


(2.6) 


(2.7) 


where  y  (=pv)  is  the  dynamic  viscosity  of  the  condensate.  Integration  of 
(2.7)  with  boundary  conditions 


w  =  0  at  x2  =  0  , 
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Here  S3  is  the  (known)  shear  stress  (in  z-direction)  on  the  condensate- vapor 
interface,  and  6(x1,z)  is  the  condensate  film  thickness.  Similarly, 
integration  of  (2.5)  with  boundary  conditions 


0  at  x. 


=  S1  at  x,  =  6 


(2.10 


yields 
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(2.11 


Here  S1  is  the  (known)  shear  stress  (in  x ,-direction)  on  the  condensate-vapor 


interface . 


It  is  advantageous  to  use  the  integral  form  of  the  continuity  equation 
rather  than  the  differential  form  (eqn.  (2.1)).  We  therefore  consider  the 
control  volume  of  Fig.  U,  and  write  the  continuity  equation  as 


A  l  *  k  / 


where  kf  is  the  thermal  conductivity  of  the  fluid,  and  hfg  is  the  latent  heat 
of  vaporization.  The  right  side  of  eqn.  (2.12)  represents  the  rate  at  which 
condensation  of  vapor  takes  place.  The  temperature  gradient  3T/3x„  can  be 


approximated  by  (Ts  -  Tw)/5,  where  Ts  is  the  saturation  temperature  at  which 

condensation  takes  place,  and  Tw  (x^.z)  is  the  wall  temperature  at  x2  =  0. 

This  approximation  basically  assumes  one-dimensional  conduction  in  the 

condensate  film  but  since  3/3x2  »  3/Sx^  or  3/3z,  and  we  are  interested  only 

in  (3T/3X.)  at  present,  it  is  Justifiable  to  make  this  approximation. 

2  6 

Substituting  for  w  and  u  from  eqns.  (2.9)  and  (2.11)  into  eqn.  (2.12), 
and  integrating,  we  get 
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The  pressure  p  can  be  related  to  the  surface  tension  o  and  radius  of  curvature 
Ri  of  the  condensate  vapor  interface  by 


p  =  ps  ±  o/Ri  , 


(2.14) 


where  p3  is  the  saturation  pressure  of  the  vapor.  The  positive  sign  holds  for 
0  <  x1  <  Xp/2,  and  the  negative  sign  for  xp/2  <  x:  <  xp,  where  xp  is  the 
length  (curve  DE  in  Fig.  3)  along  the  fluted  surface  (in  Xj-direction)  from 
the  crest  to  the  trough  (over  half  the  pitch).  We  are  basically  interested  in 
analyzing  flute  shapes  that  are  symmetric  about  the  crest  and  trough,  and 
therefore  take  advantage  of  the  symmetry  by  considering  only  half  of  the  flute 
pitch. 


Assuming  a  to  be  constant,  the  pressure  gradient  nay  be  calculated  from 


eqn.  (2.lU) 


(2.15) 


It  is  thus  the  variation  of  curvature  of  the  condensate-vapor  interface  that 
is  of  significance,  not  the  curvature  itself.  In  general,  this  can  be  quite 
different  from  d(l/R)/dx^,  a  fact  neglected  by  many  previous  analyses. 
Relations  for  finding  ifl/R^J/dx^  are  given  in  Appendix  A. 

Equation  (2.13)  is  a  partial  differential  equation  for  the  condensate 
film  thickness  6(x^,z).  It  requires  the  prior  knowledge  of  shear  stresses 
and  S3  that  come  from  vapor  dynamics  and  of  Tw  (x^z)  that  come3  from  heat 
conduction  analysis  for  the  fluted  tube  and  condensate  film.  Eqn.  (2.13) 
involves  only  first  order  derivatives  w.r.t.  z  but  fourth-order  derivatives 
w.r.t.  x1  owing  to  the  presence  of  dp/dXj^  (cf.  eqn.  (2.15)  and  (A-15)  or 
(A-18)).  For  flute  shapes  that  are  symmetric  about  the  crest  and  trough,  the 
boundary  conditions  for  the  solution  of  eqn.  (2.13)  are 


0  at  Xj  =  0,  and  at  x^ 
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for  all  z; 


(2.16a) 


5(x1,0)  =  0. 

While  the  boundary  condition  at  z  =  0  is  correct,  it  is  not  practical  to  start 
the  integration  of  eqn.  (2.13)  from  z  =  0.  Therefore,  following  Stack  and 
Merkle  [2l]  ,  we  replace  the  boundary  condition  at  z  »  0  by 


.'■-.'r.jr'ir.yv  r„ .  • .  -v-\  tv  .  .  r.  T; 


where  6Q  is  the  film  thickness  at  the  initial  station  zQ,  and  is  found  by  the 
classical  two-dimensional  Nusselt  solution 
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Here  1^00  =  (0,z0)  *  Tw  (0,0). 


(2.17) 


For  non-d imens ional i zat ion ,  we  take  the  length  L  of  the  vertical  fluted 
tube  in  the  z-direction  and  length  Xp  and  6^  as  the  characteristic  lengths  in 
the  Xj ,  and  x2-directions  respectively.  Here  <5^  is  related  to  L  by  the 
classical  Nusselt  relation  (2.17)  in  exactly  the  same  manner  as  6  is  related 
to  zQ.  Thus,  we  let 


A  -  6/6  ,  X  =  x./x  ,  Z  *  z/L  .  (2.18) 
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With  this  non-d imens ional i zat ion ,  eqn.  (2.13)  and  boundary  conditions  (2.16) 
can  be  written  as 


D  12 
U1  3X 


+ 


3E3, 
3Z  ) 


3D2, 

TT 


(2.19) 


3A  _  33A 

3X  ‘  7? 


0  at  X  *  0,1  for  all  Z, 


(2.20a) 


1/4 

A(X,Z0)  =  ZQ  , 


(2.20b) 


where 
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g(p-pv)  xp  «y 


3l£l 

%  <e-Pv)xp\ 


9 


p  -  3 

E3  “  2g(p-pv)6Y  ' 

Q  -  Ife  [rc  A/2  +  3  D3  A2/4  -  (RcA  +  I^/2D3)  In  (1+D3a/Rc)]f  , 

°3  *  5T/xP  ' 

*c  ~  *Vxp  t 

p’  !  ®  ( Tf )’  1  ' 

and 

Rn  ■  Ri/Xp 

Equation  (2.19)  is  a  highly  non-linear  partial  differential  equation, 
is  solved  numerically  by  the  finite-difference  technique  as  detailed  below, 


-  for  0  <  X  <  1/2 
+  for  1/2  <  X  <  1 


3.  FINITE-DIFFERENCE  METHOD 


Equation  (2.19)  is  parabolic  in  Z  so  that  a  forward  marching  scheme  in 
the  Z-direction  can  be  used.  Thus  the  choice  of  will  affect  only  the 
region  near  zQ ,  and  its  effect  will  die  out  an  z  increases.  This  was  indeed 
found  to  be  true.  Equation  (2.19)  as  written  is  in  conservative  form.  An 
equivalent  but  non-conservative  form  can  be  obtained  by  expanding  the 
derivatives  in  (2.19)  and  dividing  by  3 A2.  In  that  non-conservative  form,  the 
equation  is  still  parabolic  in  Z  but  the  mass  is  not  identically  conserved 
when  the  equation  is  integrated  numerically  (2ll . 

We  divide  the  interval  0  <  X  <  1  into  n  equal  parts.  Taking  6X  and  6Z  as 
the  step  sizes  in  X  and  Z-directions  respectively,  and  using  backward 
differencing  in  Z  and  mixed  differencing  in  X-direction,  we  write  the 
finite-difference  form  of  eqn.  (2.19)  as 


TszY 


i,k-l 


+  D, 


*W  * 


(1-28)Q 

Tsx 7 


i*k 


-(i-e)Q 


1-1  ,k 


6  A' 


•l*l.fc+(l-28)t  !■> -(l-9)t‘ 

rao 


i-1  ,k 


+  E_ 


i.k  - 

(fiZ) 


i,k-l 


3] 


where  the  subscripts  i  and  k  represent  the  location  in  X  and  Z  directions 
respectively  (cf.  Fig.  5),  and  0  is  a  number  that  is  selected  between  0  and  1; 
0*0  corresponds  to  backward  differencing,  0  *  1  to  forward  differencing,  and 
0  =  1/2  to  central  differencing  in  X.  Since  eqn.  (2.20b)  gives  A's  for  all  X 
(i.e.,  for  all  i)  at  ZQ  (say  k  =  1),  we  can  march  forward  in  the  Z-direction. 
In  other  words,  A's  at  (k-l)  for  all  i  are  known  in  eqn.  (3«l).  Thus  (3«l) 
leads  to  a  tridiagonal  set  of  non-linear  algebraic  equations  to  be  solved  for 
(n+l)  values  of  A’s  at  each  location  k  of  forward  march  in  the  Z-direction. 
This  non-linear  set  is  solved  by  linearization  and  successive  iteration.  Let 
us  write  eqn.  (3«l)  in  the  form 


i-1  ,k 


)" 


+ 


\m 


l+l ,k/ 


(3.2) 


where  (  )m  represents  the  values  at  the  (current)  mth  iteration,  and  the 
coefficients  A^,  %,  C^,  and  depend  upon  the  method  of  linearization. 

For  linearization  three  methods  were  tried.  In  the  first  method,  we 
simply  let 


(3.3) 


where  (  )m_1  represents  the  known  values  at  the  (previous)  (m-l)th  iteration. 
The  iteration  is  started  by  taking  A^  ^  =  Aj^.x  f°r  i*  With  this 
linearization  the  coefficients  in  eqn.  (3.2)  are 
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and  the  boundary  conditions  in  (2.20a)  have  been  taken  into  account.  This 
method  works  but  has  two  drawbacks.  Its  convergence  is  slow,  and  it  requires 
heavy  under-relaxation  for  convergence.  The  latter  implies  that  it  is 
necessary  to  under-relax  the  values  of  A  after  every  iteration  according  to 

(A)m  =  (A)0’1  +  X  j^(A)m  -  (A)®-1  j  ,  (3.5)< 

where  the  relaxation  factor  X  <  1.  In  fact,  it  was  found  that  X  decreases 
from  about  0.5  at  ZQ  to  very  small  values  as  Z  increases.  Typically  X  <  0.1 
for  Z  >  0.02. 

In  the  second  method,  we  used  Taylor  series  expansion  in  order  to 
linearize  terms  containing  A2  and  A3.  Thus  we  let 


(3.6a 


and 


) 


(3.6b 


With  this  linearization  the  coefficients  in  eqn.  (3.2)  are 

Ai  *  -  §  (1-s)[Di(°i.i,k)n"1  +  2  Ei('Wf1]> 1=1-2 . (n-1)  <3-Ta 


*  This  equation  is  to  be  treated  in  the  context  of  an  equivalent  FORTRAN 
statement . 
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3  2  3D2(5Z) 
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where  G  is  again  given  by  (3«^g)  and  the  boundary  conditions  in  (2.20a)  have 
been  accounted  for. 

It  is  of  interest  to  note  from  eqns.  (3 -T )  that  the  Taylor  series 
expansion  was  not  used  to  linearize  the  right  side  of  eqn.  (3.1)  as  well  as 
the  terms  involving  Q  on  the  left  side  of  (3.1).  While  such  a  linearization 
of  the  right  side  of  (3.1)  can  be  carried  out  without  much  difficulty,  it  was 
not  found  to  be  beneficial.  The  unnecessary  complication  was  therefore 
avoided.  The  Taylor  series  linearization  of  terms  involving  Q  in  (3.1)  was 


also  tried  but  was  found  not  to  work  at  all.  This  method  was  found  to  be  the 
best  of  the  three.  It  allows  a  larger  step  size  ( <5Z)  in  the  Z-direction,  and 


a  larger  relaxation  factor  X  than  the  first  method.  Overall,  it  is  better 
than  the  first  method  by  a  factor  of  at  least  twenty,  and  was  therefore 
adopted.  The  third  method  consisted  in  the  use  of  Hewton-Raphson  method  for 
solution  of  the  non-linear  equations.  A  general  outline  of  this  method  is 
given  in  Appendix  B.  Though  this  method  does  not  require  any  relaxation 
factor,  it  produced  absurd  results  (A's  oscilling  with  X  after  a  few  step 
sizes  in  the  Z-direction) ,  and  was  therefore  abandoned. 


4.  COMPUTATIONAL  DETAILS 


Before  the  set  of  equations  (3.2)  resulting  from  eqn.(3.l)  can  be  solved, 
we  need  to  specify  values  of  the  dimensionless  parameters  E,,  E3  and  D  2«  This 
requires  the  prior  knowledge  of  shear  stresses  S ^and  S 3>  and  the  temperature 
Tw(x1,  z)  of  the  condensate-wall  interface.  The  computer  code  does  have  the 
provision  for  specification  of  shear  stresses  S1  and  S3  but  for  the  present 
results  both  these  stresses  were  set  to  zero.  This  is  because  vapor  dynamics 
that  determines  these  stresses  is  beyond  the  scope  of  the  present  work.  The 
computer  code,  however,  does  calculate  TyCx^z)  by  considering  two-dimensional 
conduction  within  the  fluted  tube  as  well  as  the  condensate  film.  For  this, 
the  Laplace  equation 


32T  \ 
2  / 


0 


(4.1) 


is  solved  subject  to  the  boundary  conditions 

^  =  0  at  x  *  0,  P/2  where  P  =  pitch  of  the  flute, 

T  =  Ts  at  the  vapor-condensate  interface,  (4.2) 

and  T  =  Tc  at  the  coolant-tube  interface, 
where  k  is  the  thermal  conductivity,  and  Tc  is  the  coolant  temperature  assumed 
to  vary  linearly  from  the  coolant  inlet  to  exit  temperature.  The  last 
boundary  condition  in  (4.2)  needs  modification  due  to  a  film  resistance  on  the 
coolant  side  but  such  a  specification  either  requires  an  ad-hoc  value  of  the 
film  resistance  or  a  consideration  of  coolant  dynamics  that  was  also 
considered  beyond  the  scope  of  the  present  analysis.  Solution  of  eqns.  (4.1) 
and  (4.2)  was  obtained  by  a  finite  element  method  in  which  the  region  0ABC0 


(Pig. 3)  was  divided  into  several  linear  triangular  elements.  Details  for  this 
can  be  found  in  any  text  on  the  finite  element  method.  Some  care  is,  however, 
required  since  the  thermal  conductivity  of  the  fluid  in  region  ABEDA  is  vastly 
different  from  that  of  the  tube  material  in  region  DECOD  (Fig.3).  The  finite 
element  method  was  preferred  over  the  finite-difference  method  owing  to  the 
irregular  shape  of  surface  AB  (Fig.  3).  From  this  solution,  we  get  the  values 
of  Tw(x1,z)  on  the  surface  DE  (Fig.3). 

Since  the  solution  of  eqn.  (3.1)  requires  Tw(x.,z)  which  comes  from  the 
solution  of  (4.1)  for  which  knowledge  of  6-values  is  essential  (for  locating 
surface  AB  in  Fig.3),  the  two  solutions  are  coupled.  Also,  since  the  solution 
of  (3.1)  is  found  iteratively,  it  implies  that  eqn. (4.1)  should  be  solved  at 
every  iteration.  This  places  a  rather  prohibitive  demand  on  the  computer 
time.  However,  since  (6Z)  is  small,  of  the  order  of  10” 5 ,  changes  in  6  are 
small  at  every  iteration  with  the  result  that  it  is  not  really  necessary  to 
solve  equation  (4.1)  at  every  iteration.  Thus  (4.1)  was  solved  only  once  per 
step  in  the  Z-direction.  This  saves  considerable  computer  time  since  the 
finite  element  solution  of  (4.1)  requires  almost  as  much  computer  time  as 
nearly  one  hundred  iterations  for  solution  of  (3.1).  Some  sample  runs  were 
initially  made  to  confirm  that  the  error  in  6-values  due  to  this  time  saving 
feature  was  really  negligible. 

4.1  Condensate-Vapor  Interface  Profile 

As  pointed  out  in  Appendix  A,  one  method  for  finding  d(l/Ejm)/dx^  is  by 
use  of  equation  (A. 8)  which  involves  derivatives  of  f^(x),  where  f^(x) 
describes  the  condensate-vapor  interface  at  any  Z.  Several  methods  were  used 
to  find  these  derivatives  but  none  of  them  proved  worthwhile  since  f^(x)  is 
known  at  unequally-spaced  values  of  x.  The  methods  are  listed  below  for  the 


sake  of  completeness. 


i)  Cubic  splines  were  fitted  in  a  least  squares  manner  through  fj_  values 
to  find  fj_'  and  f^".  Then  cubic  splines  were  fitted  through  the  previously- 
computed  f  values  to  find  f^' ' ' .  A  slight  variation  of  this  method  was  to 
fit  cubic  splines  repeatedly  through  f  ^ ,  f ^ '  and  f values  to  find  ' ' ' . 

Both  these  methods  gave  different  values  for  f^"  and  f^' ' '  pointing  to  their 
futility. 

ii)  Parabolas  were  fitted  through  three  consecutive  points  in  succession 
on  the  condensate-vapor  interface,  and  Rj.  and  d( l/R-^/dx x  were  computed  at  the 
point  central  to  each  set  of  three  points.  This  also  led  to  inconsistent 
results . 

iii)  An  attempt  to  fit  a  truncated  Fourier  series  (having  3-U  terms) 

through  f^  values  was  made.  This  was  attempted  since  we  were  working  with  a 

sinusoidal  flute.  However,  since  the  fit  itself  was  poor,  no  attempt  was  made 

to  calculate  the  derivatives.  A  better  fit  (Fig.  6)  was  obtained  by  fitting 

« 

(ax+  &2  cos  2-irx/P)  through  all  (n+l)  values  of  f*,  where  the  coefficients  ax 
and  a2  were  found  for  a  least  squares  fit.  However,  even  this  fit  was 

considered  inadequate  to  find  d(l/Ri)/dxx.  A  slight  variation  of  this  was  to 

fit  (ax+  a2  cos  2irx/P)  through  the  first  n/2  values  of  f*  (for  x  <  P/2)  and  a 
similar  curve  through  the  last  n/2  values.  This  produced  a  very  good  fit  on 
visual  inspection  (Fig.7).  However,  using  it  to  find  d(l/Ri)/dx  ^  leads  to 
wiggles  in  6-values.  Solution  of  equations  (3.1)  and  (U.l)  starting  from  ZQ 
leads  to  film  shape  rising  suddenly  in  the  trough  region  (Fig. 8)  after  a  few 
(6Z)  steps,  thus  violating  the  boundary  conditions  there.  Starting  from 
Z  *  0.07,  with  results  at  Z  *  0.07  found  by  better  methods,  the  solution  leads 
to  such  large  wiggles  that  even  negative  values  of  6  are  found  only  a  few 

steps  downstream.  In  an  effort  to  improve  upon  this  situation,  values  of  5 

were  computed  to  correspond  to  the  fit  in  Fig.  7  before  iterating.  This, 


however,  slowed  down  the  convergence  so  much  that  even  two  hundred  iterations 
were  insufficient  for  the  same  tolerance.  In  fact,  it  appeared  as  if 
convergence  would  never  take  place.  Nevertheless,  this  exercise  did  point  out 
the  reason  for  the  failure  of  the  above  technique,  since  the  differences 
between  the  original  5-values  and  those  conforming  to  the  good  ( ! )  fit  in 
Fig.  7  were  found  to  be  large,  as  much  as  20%. 

iv)  Efforts  to  use  lagrangian  interpolation  to  compute  f^  values  at 
equi-distant  x  values,  and  then  use  techniques  described  in  Section  h.2  to 
find  f^' ,  etc.  also  produced  absurd  results.  Use  of  equation  (A. 8)  to 
compute  d(l/Ri)/dx1  was  therefore  abandoned. 

It  may  be  mentioned  that  while  using  eqn.  (1.8)  to  find 
dfl/R-j^/dXj ,  a  futile  attempt  was  also  made  to  solve  the  partial  differential 
equation  (2.13)  directly  in  terms  of  63  rather  than  6. 

U.2  Derivatives  of  Condensate  Film  Thickness 

Use  of  either  eqn.  (A. 15)  or  (A.l8)  to  find  d(l/Ri)/dx1  requires  the 
determination  of  first  three  derivatives  of  5  w.r.t.  x^.  Fortunately  5  is 
known  at  equidistant  values  of  x^.  It  is  therefore  easier  to  find  these  der¬ 
ivatives  than  those  of  f^  w.r.t.  x.  Several  methods  for  this  determination 
were  also  tried  with  varying  degrees  of  success. 

i)  Derivatives  of  delta  were  found  by  fitting  a  cosine  to  the  vapor- 
condensate  interface  in  a  manner  analogous  to  that  described  in  Sec.  U.l,  but 
it  resulted  in  absolutely  useless  results. 

ii)  Another  fruitless  effort  was  to  fit  least  squares  cubic  splines, 
repeatedly  or  otherwise,  through  5-values  so  as  to  find  s' ,  5"  and  s' ' ' . 


Though  these  sets  of  linear  equations  are  tridiagonal  and  can  be  solved 
easily,  they  do  represent  an  increase  in  computer  time.  Besides  the  method 
yields  rather  poor  results. 


iv)  A  method  that  gives  better  results  than  those  obtained  by  above 
methods  is  based  on  passing  least  squares  quartics  through  seven  consecutive 
points  at  a  time.  Leaving  the  details  to  the  reader  (see  also  (2U]  ,  p.l*92), 

t  i> 

we  present  the  equations  for  finding  6^,  <5^  and  <5^  for  i  =  o,  1,...,  n. 

These  equations  make  use  of  the  boundary  conditions  (2.16a),  and  the  fact  that 
5  is  known  at  equidistant  x^  values. 


Relations  for  5',  accurate  to^O^dx  are 

<  '  67ai-2-  58si-i+  58  W  67  W  22 si.3)'  252  N  • 


i  =  3 *•••»( n-3 ) , 


d'  »  0  ,  6'  -  0, 

o  *  n  * 


61  =  ("586o  '  6761+  80 62  +  67 63  ~  22&k)/  252 6xi  ’ 

=  (-67dQ  -  36d1+  58«3  +  676u  -  22^)/  252dX;L  , 

*n-2  "  -  6lSn-U  *  585„-3  +  36tn-l  *  S7S„V  2525lll  - 

and 

i  “  (225  ,  -  67d  _ ,  -  80fi  ♦  67  d  ,  +  58  d  )/  2526x,  . 


Relations  for  5",  accurate  to  0 


'•Xv 


Use  of  these  equations  for  finding  5’,  6"  and  5"*  leads  to  mild 

fluctuations  in  6-values  after  a  few  steps  in  the  Z-direction.  Their  use  had 

therefore,  to  be  abandoned.  Less  accurate  relations  involving  an  error  of 

0  6^  also  yielded  similar  results. 

1 

v)  The  best  method  for  finding  the  derivatives  of  6  for  this  problem 
turned  out  to  be  the  use  of  central  difference  formulae 


f 


2  6xx 


i  =  1,2 , . . . ,(n-l) 


(4.8 


for  the  first  derivative,  and 

51  *  (  S1H  ‘  2  Si  *  Sl-l)  /  5xl  •  1-1 ’2 . (n-1)  • 

5o  ■  ( *1  -  «o  )  /  2Sxl  •  "‘•9 

d"  /«.-«)/  2  6x? 

n  *  y  n-1  n  /  1  * 


for  the  second  derivative.  The  third  derivative  was  found  by  applying  (4.8) 
to  6"  calculated  from  (4.9).  This  method  gave  the  best  results  for 
6-distributions,  and  was  therefore  adopted. 


5.  RESULTS  AND  DISCUSSION 


The  complete  computer  program  was  developed  in  two  stages.  In  the  first 
stage  the  solution  of  eqn.  (3.1)  alone  was  attempted.  This  program  was  tested 
for  the  vapor  and  tube  properties  for  which  results  are  available  in  [2l] .  As 
pointed  out  in  "Introduction",  Stack  and  Merkle  solve  a  much  simplified  form 
of  (3.1)  since  besides  other  simplifications,  they  use  Cartesian  coordinates 
rather  than  curvilinear  coordinates .  In  fact ,  we  had  modified  our  earlier 
program  to  solve  Stack  and  Merkle' s  equation  (l4)  for  the  condensate  film 
thickness,  and  obtained  identical  results.  The  finite  element  program  for  the 
solution  of  eqn.  (4.1)  was  then  developed, ’and  the  two  were  coupled. 

In  the  following  we  present  results  for  one  vapor-tube  combination  for 
which  experimental  data  is  also  available  (10,11]  for  comparison.  Since  our 
analysis  assumes  the  radius  of  curvature  of  the  fluted  tube  to  be  large 
compared  to  the  film  thickness,  it  is  not  realistic  to  consider  sinusoidal 
flutes  with  amplitude-to-pitch  ratio  much  above  0.2.  With  this  in  mind,  we 
computed  results  for  tube  'F'  of  Ref.  I 10,11].  This  tube  has  seven  skirts  and 
an  I.D.  *  22.9  nm.  Details  of  the  actual  shape  of  the  flute,  and  its 
dimensions  are  not  available  in  1 10 ,11 ] ,  and  were  therefore  approximated  from 
a  blow-up  of  the  tube  'F*  photograph  in  111] .  The  relevant  fluid  (R-113)  and 
tube  (aluminum)  properties  are 


p  =  1498.343  kg/m3 
v  =  3 .2067x10“ 7  m2/s 
hfg  *  145225.56  J/kg 
kf  *  0.06951  W/m-K 
^woo  =  318.5  K 


Py  =  8.58628  kg/m3  , 

0  =  0.01432  N/m  , 

Ts  =  325.5  K 
kv  =  205  W/m-K 

Tc  =  318.5  K  (in),  318.8  K  (out)  , 
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L  »  0.142875  m  ^  x  f  ,  ht  -  0.88077  ran  , 

P  =  1.6l4l  nsn  ,  hj  a  0.15543  mm  , 

where  kf  and  ky  are  the  thermal  conductivities  of  the  condensate  and  fluted 
tube  respectively,  and  h-t  and  ho  are  associated  with  the  flute  shape  (Fig.  3) 
taken  to  be 

f(x)  «  h^  +  hQ  cos  .  (5.1) 

Thus  ho  is  half  the  amplitude  of  the  flute.  The  fluid  properties  were 
calculated  at  322  K  (mean  of  Ts  and  Twoo)  from  known  correlations  [25].  These 
values  lead  to 

6y  -  0.08051  ran  ,  xp  =  0.876487  mm 

After  some  numerical  experimentation  involving  different  step  sizes, 
etc.,  the  solution  of  eqn.  (3«l)  was  started  from  ZQ  =  5x10“ 6  with  0  *  1/2 
(corresponding  to  central  differencing  in  X),  5X  =  0.05,  an  initial  fiZ  * 
5xl0"6,  and  an  initial  X  *  0.5  (see  eqn.  (3*5) )•  Further  iteration  for  the 
solution  of  (3.1)  was  terminated  when 

|l  "  |  kj10"1  /(Ai  kjm  |  <  £  for  411  1  at  every  k  »  (5*2) 

where  e  was  taken  to  be  10"6.  As  we  marched  downstream  in  the  Z-direction ,  <SZ 
was  increased  and  the  relaxation  factor  X  had  to  be  decreased  according  to 
Table  1,  which  also  gives  an  idea  of  the  number  of  iterations  before  (5 .2)  was 


satisfied.  Clearly,  it  takes  longer  to  get  the  converged  solution  as  we 
proceed  downstream  in  the  Z-direction.  Reasons  for  this  will  be  apparent 
soon. 


6Z 

Table  1 

Some  Parameters 

for  Solution 

of 

(3.1) 

X 

Z 

# 

of  iterations 

5xlO"6 

0.5 

up  to 

0.002 

~  15 

8x10- 6 

0.2 

up  to 

0.01 

~  30 

1.25x10“ 5 

0.1 

up  to 

0.0225 

~  50 

2xl0-5 

0.06 

up  to 

0.086 

~  90 

Three  divisions  of  the  region  OABCO  (Fig.  3)  into  finite  elements  for  the 
solution  of  eqn.  (4.1)  were  tried.  The  one  selected  on  the  basis  of  adequate 
computational  accuracy  and  relatively  economical  calculation  had  154 
triangular  elements  with  108  nodes.  There  were  11  nodes  each  on  the  faces  AB, 
BC,  OC  and  DE  (Fig.  3),  and  15  nodes  on  the  face  AO. 

Fig.  9  shows,  to  true  scale,  the  flute  shape  (thick  curve),  and  the 
condensate  film  shape  (dashed  curves)  at  four  values  of  Z  =  0.001,  0.01,  0.03, 
and  0.086.  As  expected  the  film  thickens  quite  rapidly  in  the  trough  region 
while  it  remains  thin  over  the  crest.  Unfortunately  we  are  unable  to  present 
any  results  for  Z  much  greater  than  0.086  since  convergence  of  our  solution 
for  the  film  thickness  (with  the  same  tolerance  e)  becomes  very  slow.  The 
principal  reason  for  this  slow  convergence  is  the  thickening  of  the  film  in 
the  trough.  Moreover,  it  is  really  futile  to  try  to  speed-up  the  convergence 


since  our  formulation  for  the  condensate  film  thickness  is  on  shaky  ground  as 
the  film  thickens.  The  non-linear  partial  differential  equation  (2.13)  comes 
from  a  boundary-layer  type  analysis ,  i.e.,  we  neglect  all  velocity  gradients 
except  those  in  the  direction  normal  to  the  fluted  surface.  While  this  is 
perfectly  reasonable  over  the  initial  portion  of  the  tube  height  when  the 
condensate  film  is  thin,  it  is  increasingly  questionable  as  the  film  thickens 
in  the  trough  region.  It  appears  that  for  the  fluid-tube  combination 
considered  here,  it  is  improper  to  analyze  the  film  in  the  trough  region  in 
the  boundary-layer-sense  for  Z  >  0.085*  It  is  also  clear  from  Fig.  10,  that 
shows  the  absolute  value  of  d(l/R^)/dxi  at  various  Z  values,  that  this 
derivative  is  negligible  in  the  trough  region  as  compared  to  it3  -value  in  the 
crest  region.  It  is  therefore  appropriate  to  neglect  the  surface  tension 
effects  in  the  trough  region  once  the  condensate  film  has  thickened.  Fig.  10 
also  plots  the  absolute  value  of  the  derivative  of  flute  curvature  with  x^, 
and  even  at  the  low  value  of  Z  *  0.002,  it  is  very  different  from  the 
derivative  of  condensate-vapor-interface-curvature  with  X]_. 

While  these  considerations  will  be  undertaken  in  an  extension  of  this 
project,  it  is  encouraging  to  note  from  Fig.  11  that  a  linear  extrapolation, 
on  a  log-log  basis,  of  our  results  to  date  shows  a  fairly  good  comparison  with 
experimental  data.  Fig.  11  shows  the  heat  load  per  half  flute  (in  W)  as  a 
function  of  Z  on  a  log-log  scale.  The  solid  portion  of  the  straight  line  is 
based  on  computed  results  while  the  dashed  portion  is  the  extrapolation.  The 
lone  circled  point  is  based  on  the  experimental  data  [ 11 ] .  We  also  checked 
that  heat  transfer  rates  across  the  faces  AB  and  0C  (Fig.  3)  matched  within 

Theoretically  they  should  be  identical  since  faces  0A  end  BC  are 


insulated 


(I 
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6.  SUGGESTIONS  FOR  FURTHER  WORK 

Besides  the  deficiencies  mentioned  above,  our  present  analysis  also 
neglects  the  transverse  curvature  term,  i.e.,  it  assumes  the  film  thickness  to 
be  small  compared  to  the  radius  of  curvature  of  the  flute.  The  analysis  is 
therefore  restricted,  in  its  present  form,  to  flutes  having  a  small  amplitude 
to  pitch  ratio.  This  is,  however,  a  minor  drawback  since  it  seems 
(analytically  at  least)  that  we  can  overcome  this  deficiency  relatively 
easily.  Its  implementation  in  the  computer  code  is  yet  to  be  done. 

Another  assumption,  common  to  all  analyses  to  date,  and  far  more  serious 
in  terms  of  correct  modelling  of  the  practical  applications,  is  the  complete 
neglect  of  vapor  shear  on  the  interface.  Our  code  does  have  the  provision  for 
studying  its  effect  but  since  we  have  not  analyzed  the  vapor  dynamics  yet,  we 
simply  take  the  vapor  shear  to  be  zero  at  present.  One  should  also  consider 
the  coolant  dynamics  for  a  complete  solution  but  all  this  will  undoubtedly  be 


very  demanding 
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APPENDIX  A 


CURVATURE  OF  CONDENSATE-VAPOR  INTERFACE 
AND  ITS  VARIATION  ALONG  FLUTE  SHAPE 

Equation  (2.13)  for  the  determination  of  condensate  film  thickness 

6(x  ,z)  involves  the  value  of  pressure  gradient  dp/dx^  which,  in  turn, 

depends  upon  d(l/R^)/dx^,  as  is  clear  from  (2.15).  We  present  below  some 

relations  for  calculation  of  this  quantity. 

The  basic  terminology  is  present  in  Fig.  3,  wherein,  f(x)  denotes  the 

flute  shape  and  f^(x)  the  condensate-vapor  interface.  The  position  vector  to 

a  point  on  the  flute  is  rw,  and  that  to  a  point  on  the  condensate- vapor 
-► 

interface  is  r^,  such  that 

ri  -  rw  =  ny.5  ,  (A.l) 

where  n^.  is  the  unit  normal  to  the  flute  as  shown  in  Fig.  3.  We  also  let  x^ 
be  the  curvilinear  coordinate  along  the  condensate-vapor  interface. 

Then  by  definition 


1  _  a2ri 

Ri  “  axj2 

(A.2) 

Now  rj^  *  (x,  fi). 

(A.3) 

x  »? 

and  xt  *  £ (l  +  f*  )  d?  , 

(A.U) 

where  prime  denotes  differentiation  with  x.  From  equation  (A.U),  we  get 


While  equation  (A. 8)  seems  to  be  in  a  convenient  form  for  calculation  of 


d(l/Ri)/dx^,  it  is  not  practical  since  f^  is  known  only  at  a  discrete  set  of 
unequally-spaced  x-values.  Thus,  finding  fj.'  ,  f^"  and  fj_" '  is  troublesome,  to 
say  the  least.  Several  attempts  to  use  equation  (A. 8)  had  to  be  abandoned  for 
this  reason  (see  Sec.  4.1  for  details). 

Another  relation  for  the  curvature  of  the  condensate-vapor  interface 
comes  from  use  of  equation  (A.l),  and  the  fact  that  the  unit  normal  vector 
is  given  by 


(  dx  % 

dxx  *  dxf  1 


■  (1+f 
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-lfZ 
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(-f',  1)  . 


(A. 10) 


Thus,  from  (A.l) 
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(1+f'2) 


f  + 


(1+f'2) 


V2 


j  ,  since  rw  =  (x,f ). 


dr 


Therefore , 


dx 


-  *  ( i  -  «v. 


6f 


(If*2) 
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,  t'+  $'  - 


6f '  f  " 


(1+f  ) 


3/z 


(A. 11) 


where  6  =  36/3x^  but  f'=  df/dx,  etc..  The  reason  for  keeping  36/3x^  in  place 


of  36/3x  is  that  solution  of  equation  (2.13)  leads  to  values  of  6  at  equidis¬ 


tant  x^-values  but  at  non-equidistant  x-values.  It  is  therefore  convenient 


to  find  numerically  derivatives  of  6  with  respect  to  x  . 


It  is  interesting  to  note  that  many  theoretical  analyses  that  claim  to 
determine  the  surface  tension  effect  based  upon  the  actual  shape  of  the 
condensate-vapor  interface  invariably  use 


*i  (1+6'V/2 


(A.lU) 


for  the  curvature  of  the  interface.  This  follows  from  the  correct  relation 
(A. 13)  if  only  the  first  terms  in  both  the  numerator  and  the  denominator  are 
retained.  However,  even  for  small  d,  there  is  hardly  any  justification  for 
dropping  all  other  terms  in  (A.13).  The  difference  between  the  actual 
d(l/Rj[)/dx^  and  that  obtained  from  (A.lU)  gets  further  amplified  due  to 
differentiation. 

For  convenience  let  us  write 


1 

*i 


N.s 
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gn(N) 

V2 
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where  the  expressions  for  N  and  D  are  clear  from  equation  (A.13).  It  is  then 
easy  to  get 


sgn(W) 

i/z 
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(V5V 


+  «  (T= 


«<vV> 


}  +  Yi 


+  6  2(y,i+  2y1'  )  +  6'(y3+  5(2y2+  Y^')}  +  52y2’  +  «Y3' 


-  2rJd“  ^6'(2^66"  +  P'jS  +  V  +  (1  +  5 * 2 ) T6 '  +  S2Y7'  +  «Y8'},  (A. 15) 


Equation  (A. 15)  is  quite  complicated  but  as  mentioned  in  the  discussion  on 

3  3 

eqn.  (2.13),  it  contains  3  6/3x^  .  It  is  therefore  essential  to  determine  the 
condensate  film  thickness  as  accurately  as  possible.  Otherwise,  errors  in 
finding  the  third  derivative  of  6  w.y.t.  x^  numerically  will  be  excessive. 

Due  to  these  numerical  difficulties,  it  was  found  that  use  of  eqn.  (A. 15)  to 
find  d(l/R^)/dxi  resulted  in  mild  oscillations  in  6-values  after  a  few  steps 
in  the  Z-direction.  Several  methods  (at  least  five)  were  used  to  find 
6-derivatives  w.r.t.  x^  but  all  resulted  in  oscillatory  6-values;  the  least 
oscillations  resulted  from  the  use  of  central-difference  formulae. 

In  order  to  overcome  this  difficulty,  another  relation  for  1/R^  and  its 
derivative  w.r.t.  x^  was  developed.  This  relation,  based  on  eqns.  (A. 2)  and 
(A.l),  is 


32(rw  +  n^S) 


32?w  32^ 
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3xi2 


3^w 
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36  •*  32  6 

3x.  "w  2 

1  Sx^ 


(A.16) 


where  3/3x^  has  been  approximated  by  3/3x^.  Noting  that  rw  =  (x,f),  and  nw 
is  given  by  eqn.  (A.10),  we  get  from  (A.l6),  after  some  algebra,  that 
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where  6  =  35/3x^  ,  f  =  df/dx,  etc.,  and 
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From  equation  (A. 17)  it  follows  that 
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Equation  (A.l8)  also  contains  3  6/3x^  ,  but  unlike  equation  (A. 15),  it  yields 
stable  (non-oscillatory)  values  of  6  upon  integration  of  eqn.(2.13).  This 


relation  was  therefore  used  in  the  final  computer  code. 


APPENDIX  B 


Nonlinear  Equations  Solution  by  Nevton-Raphson  Method 


Consider  the  set  of  equations 


(K!  {x}  *  {R} 
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2  kij(x,  ’  V  **’  *n^xJ  =  ^i^x-,  ’  V  **’  1 

!=1  1  2  12 


1 ,  ■  n a 


where  k^j  and  R^  are  non-linear  functions  of  x  ,  x  ,  Xp. 

n 

Let  F^Cx^,  x^,..,  x,^)  =  I  k^j(x^,  x^,  -  R^(x^,  x^,  xn), 

X  =  1)  2 ,  a  a  ,  Ha 


So  is  the  residue  in  the  i  eqn.  For  solution  Fi  =  0,  i  =  1,  2,  n. 
Expanding  F^  by  Taylor  series,  we  write 


F.i  ( x  +  Ax  ,  x  +  Ax  ,  . . ,  x_  +  Axn ) 
1  12  2  n 
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12  J»1  '5x7  1  2  n  J 


i  =  2 )  •  • ,  n. 


Neglecting  all  higher-order  terms  and  setting  leftside  =  0  for  a  solution, 
we  get 
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Thus  Newton-Raphson  iteration  algorithm  can  be  written  as 


,  ,  m  m+1  m 

l J]  {Ax}  =  -{F> 


m+1  m  m+1 

{x}  =  {x}  +  {Ax} 


.  m 


where  m  denotes  mth  iteration  &  (j)  is  the  system  Jacobian  at  the  m™  itera¬ 


tion. 


m 


The  elements  of  |J]  are 


Jij  -  ^7  (xT  »  X2“  ,  ••  .  xR) 
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If  it  is  too  time  consuming  to  compute  (j]  at  every  iteration,  find  [j]  once 
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(b)  Cross  section  of  vertical  condenser. 
Fig.  1.  Vertical  Submarine  Condenser 
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